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Abstract — The accuracy of different transfer matrix ap- 
proaches, widely used to solve the stationary effective mass 
Schrodinger equation for arbitrary one-dimensional potentials, 
is investigated analytically and numerically. Both the case of a 
constant and a position dependent effective mass are considered. 
Comparisons with a finite difference method are also performed. 
Based on analytical model potentials as well as self-consistent 
SchrSdinger-Poisson simulations of a heterostructure device, it 
is shown that a symmetrized transfer matrix approach yields a 
similar accuracy as the Airy function method at a significantly 
reduced numerical cost, moreover avoiding the numerical prob- 
lems associated with Airy functions. 

Index Terms — Quantum effect semiconductor devices, Quan- 
tum well devices, Quantum theory, Semiconductor heterojunc- 
tions, Eigenvalues and eigenfunctions, Numerical analysis, Tun- 
neling, MOS devices 

I. Introduction 

TRANSFER matrix methods provide an important tool 
for investigating bound and scattering states in quantum 
structures. They are mainly used to solve the one-dimensional 
Schrodinger or effective mass equation, e.g., to obtain the 
quantized eigenenergies in quantum well heterostructures and 
metal-oxide-semiconductor structures or the transmission co- 
efficient of potential barriers (TJ, 0, 0, 0. Analytical 
expressions for the transfer matrices are only available in 
certain cases, as for constant or linear potential sections and 
potential steps 0. An arbitrary potential can then be treated 
by approximating it for example in terms of piecewise constant 
or linear segments, for which analytical transfer matrices 
exist. For constant potential segments, the matrices are based 
on complex exponentials (TJ, 0, while the linear potential 
approximation requires the evaluation of Airy functions 0. 

Many applications call for highly accurate methods, e.g., 
quantum cascade laser structures where layer thickness 
changes by a few A already lead to significantly modified 
wavefunctions, resulting in altered device properties 0, |6). 
Also numerical efficiency is crucial, especially in cases where 
the Schrodinger equation has to be solved repeatedly. Exam- 
ples are the shooting method where the eigenenergies of bound 
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states are found by energy scans, or Schrodinger-Poisson 
solvers working in an iterative manner [3]. Besides providing 
accurate results at moderate computational cost, an algorithm 
is expected to be numerically robust, and a straightforward 
implementation is also advantageous. 

Besides transfer matrices, also other methods are frequently 
used, in particular finite difference or finite element schemes 
0, 0. For scattering state calculations, they are compli- 
mented by suitable transparent boundary conditions, resulting 
in the Quantum Transmitting Boundary Method (QTBM) 0, 
|9l . The transfer matrix method tends to be less numerically 
stable than the QTBM, since for multiple or extended barriers, 
numerical instabilities can arise due to an exponential blowup 
caused by roundoff errors 0. This issue can however be 
overcome, for example by using a somewhat modified matrix 
approach, the scattering matrix method iflOll . In this case, the 
transfer matrices of the individual segments are not used to 
compute the overall ttansfer matrix, but rather the scattering 
matrix of the structure. In addition, transfer matrices have 
many practical properties, such as their intuitiveness partic- 
ularly for scattering states, the intrinsic current conservation, 
and the exact treatment of potential steps, which arise at the 
interfaces of differing materials. This makes them especially 
suitable and popular for 1-D heterostructures or metal-oxide- 
semiconductor structures, providing a simple, accurate and 
efficient simulation method 0. 

As mentioned above, transfer matrices are usually based 
on a piecewise constant or piecewise linear approximation 
of an arbitrary potential, giving rise to exponential and Airy 
function solutions, respectively. The main strength of the Airy 
function approach is that it provides an exact solution for 
structures consisting of piecewise linear potentials, and hence 
only requires few segments for approximating almost linear 
potentials with sufficient accuracy. On the other hand, Airy 
functions are much more computationally demanding than 
exponentials, and also prone to numerical overflow for regions 
with nearly flat potential QT|. Thus, great care has to be taken 
to avoid these problems, and to evaluate the Airy functions in 
an efficient way fl2l . 

It would be desirable to combine the advantages of both 
methods, namely the accuracy of the piecewise linear approx- 
imation and the computational convenience of the exponential 
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Fig. 1. Various transfer matrix schemes applied to segmented potential. 
Shown is the exact (solid line) and approximated (dashed line) potential, (a) 
Piecewise constant potential approximation, (b) Piecewise linear approxima- 
tion, (c) Piecewise contant approximation for symmetrized transfer matrix. 



transfer matrix scheme. In this paper, we evaluate the accuracy 
and efficiency of the different transfer matrix approaches, 
taking into account both bound and scattering states. In this 
context, analytical expressions for the corresponding local 
discretization error are derived. We furthermore evaluate the 
different approaches numerically on the basis of an analytically 
solvable model potential, and also draw comparisons to the 
QTBM. In particular, we demonstrate that a symmetrized 
exponential matrix approach is able to provide an accuracy 
comparable to that of the Airy function method, without 
having its problems and drawbacks. In our investigation, we 
will consider both the case of a constant effective mass and 
the more general case of a position dependent effective mass. 

II. Transfer matrix approach 

In a single-band approximation, the wavefunction ip of 
an electron with energy £ in a one-dimensional quantum 
structure can be described by the effective mass equation 

h2 -d z —^d z + V{z)-E\ i){z) = 0. (1) 



2 ~"m* (z) 

Here, the effective mass m* and the potential V generally 
depend on the position z in the structure. For applying the 
transfer matrix scheme, we divide the structure into segments, 
see Fig. [T] which can vary in length. Potential and effective 
mass discontinuities can be treated exactly in transfer matrix 
approaches by applying corresponding matching conditions. 
To take advantage of this fact and obtain optimum accuracy, 
the segments should be chosen so that band edge discon- 
tinuities, as introduced by heterostructure interfaces, do not 
lie within a segment, but rather at the border between two 
segments. 

A. Conventional transfer matrices 

For the piecewise constant potential approach (Fig. [Ha)), the 
potential and effective mass in each segment j are approxi- 
mated by constant values, e.g., Vj = V (Zj), m* = m* (zj) 
for zj < z < zj + Aj 



Zj+i, and a jump Vj — > Vj+%, 



([T|l is for z-j < z < z 



at the end of the segment J2j. The solution of 



3+1 



then given by 



ip (z) = Aj exp [ikj (z — Zj)] + Bj exp [— ikj (z — Zj 



(2) 



where kj = j2ra*(E—Vj)jh is the wavenumber (for 



E < Vj, we obtain % = iKj = iW2mj (Vj - E)/H) 0. 
The matching conditions for the wavefunction at the potential 
step read 

4> (ZQ+) = "0 (Z0~) , 

[d z i> ( 2o +)] /m* (z +) = [d z i> (z -)] /m* (z -) , (3) 

where zq+ and zq— denote the positions directly to the right 
and left of the step, here located at zq = Zj+\ 0. The 
amplitudes Aj + i and -Bj+i are related to Aj and Bj by 



3 + 1 
3 + 1 



T 



3,3+1 



Aj 
Bi 



(4) 



with the transfer matrix 



Tj,j+1 — T i^j+~L T 3 (A?) 



%+i+fe r ifc,A,- 
2ft + l 



3 J + i+ft- r -ifc,A, ' ( 5 ) 
2ft + i C / 



Equation (0 is the product of the transfer matrix for a flat 
potential 







obtained from (O, and the potential step matrix 

03+1+ 0j 03 + 1-03 



T, 



1 



3->3+l 



2ft 



3+1 



0. 



3 + 1 



03 03 +1 +03 



(6) 



(7) 



with ftj — kj/m*, derived from (0! J4). The relation between 
the amplitudes at the left and right boundaries of the structure, 
Aq, Bo and An, Bn, can be obtained from 



Bn 



N-1,nTn-2,N-1 



■ T, 



0,1 



A 
Bo 



Tn T12 
T21 T22 



Ao 
Bo 



(8) 



where N is the total number of segments. For bound states, 
this equation must be complemented by suitable boundary 
conditions. One possibility is to enforce decaying solutions at 
the boundaries, Aq = Bn = 0, corresponding to T22 = 
in (0, which is satisfied only for specific energies E, the 
eigenenergies of the bound states J2|- 

For the piecewise linear potential approach (Fig. 1(b)), the 
potential in each segment j is linearly interpolated, V (z) = 
Vj + V z .j (z — Zj) for Zj < z < Zj + Aj = Zj+i, with V z .j — 
(Vj+i — Vj) / Aj. Equation ((TJ can then be solved analytically 
in terms of the Airy functions Ai and Bi (2), 



i>(z) = Aj Ail Sj + 



BjBi 1 *, 



£3 



(9) 
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for Zj < z < Zj+i, with Sj = (Vj — E) /ej and lj = £j/V Zl j, 
where Sj = ytPV^j] (2m*). We obtain 



and 



A, 



A, 



(10) 



(11) 



'3 *"\°3Jr3 1 "3 *'J^ 1 V3/rj, 

with Dj = Ai(sj)Bi'(sj) — Ai' (sj)Bi(sj). Here a prime 
denotes a derivative with respect to the argument of the Airy 
function (for Ai', Bi') or the position z (in all other cases). 
A position dependent effective mass is treated by assigning 
a constant value to each segment j, for example m* (zj) 
or preferably [to* (zj) + to* (zj + i)] /2 (see appendix), and 
using the matching conditions Q at the boundary between 
two adjacent segments (2). A piecewise linear interpolation 
of to* as for the potential is not feasible, since then the 
solutions of (Q~|i cannot be expressed in terms of Airy functions 
anymore. Equations ( [Tol l, (fTTT i can again be rewritten as a 
matrix equation of the form ©, allowing us to treat the 
quantum structure using (© in a similar manner as described 
above [2|. Interfaces introducing abrupt potential changes in 
the quantum structure must be taken into account explicitly 
in the Airy function approach by employing the matching 
conditions (0. 

B. Symmetrized matrix 

In the transfer matrix approach, the amplitudes Am and 
at the right boundary of the structure are related to the values 
A and B at the left boundary by repeatedly applying the 
transfer matrix. Due to the segmentation of the potential, an 
error is introduced in © for every propagation step from a 
position zj to Zj+i, which is typically characterized in terms 
of the local discretization error (LDE). The LDE is defined as 
the difference between the exact and computed solution at a 
position Zj + i obtained from a given function value at Zj. In 
the appendix, the LDE with respect to the amplitudes Aj and 
Bj for the transfer matrix <(5j is found to be O (A|). It can be 
improved to O (A|) by symmetrizing the matrix, i.e., placing 
the potential step in the middle of the segment, see Fig. 1(c). 
The resulting transfer matrix is then with — {kj ± kj+i) /2 
given by 



T 3,j+1 — T 3 + ^ 



2) T ^+ lT * IT 



ft+i+fty ifc+A,- 

2ft + l e 



ft+i+fey -ifc+A,- 
2ft + i L 



(12) 



As in the 



where again kj — \j2rn* (E — Vj)/h, f3j — kj/m*. 
Airy function approach, interfaces introducing abrupt potential 
changes in the quantum structure must be dealt with separately 
by applying the matching conditions; here, the corresponding 
transfer matrix (|7]) can be used. 
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Fig. 2. Exponential model potential with d = 10 nm and K = — l/d, used 
for evaluating the accuracy of various methods, (a) Barrier, (b) Quantum well. 



III. Comparison 

The improved transfer matrix ( fT2l can be evaluated at a 
comparable computational cost as the matrix (0, but exhibits 
a superior accuracy. As shown in the appendix, the local 
discretization error with respect to the amplitudes Aj and Bj 
is improved from O (A|) to O (Aj) for arbitrary potentials 
and effective masses, i.e., the same order as for the Airy 
function approach, which however involves a significantly 
higher computational effort. 

In the following, we compare the accuracy of the different 
methods for an analytically solvable model potential. Here, 
polynomial test potentials are not suitable for a general dis- 
cussion since their higher order derivatives identically vanish, 
which can lead to an increased accuracy in such special cases. 
Especially triangular or other piecewise linear potentials are 
obviously inadequate since the Airy function approach then 
becomes exact. Instead, we choose the exponential ansatz 

V(z) =V + Vi exp (Kz) , (13) 

< z < d (see Fig. |2), approaching a linear function for 
K — > 0. Such a potential can for example serve as a model for 
the effective potential profile in the presence of space charges 

ma, d. 



A. Position independent effective mass 

For now, we assume a constant effective mass to* 
analytical solutions of the form 



Then, 



(14) 



ip = ci J M (a) + c 2 Y ll (a) 

exist for the potential ( fT3l , with constants c\ and C2. Here, J M 
and Y M are Bessel functions of the first and second kind, and 
the parameters are given by 



/' 

a(z) 



, V2to* (Vq - E) 
KK 



1 



(15) 



For our simulations, the different transfer matrix approaches 
discussed in Section [TT] are used to compute an overall matrix 
based on ([8]), from which the required quantities can be 
extracted. First, we investigate the barrier structure shown in 
Fig. |2 a), which can be characterized in terms of a transmission 
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Fig. 3. Relative error et = |1 — T nuul /T\ of the numerically obtained 
transmission coefficient T num as a function of the number of segments TV. 
The corresponding barrier is shown in Fig. 12a), the effective mass is assumed 
to be constant. 
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Fig. 4. Relative error et = |1 — T n um/T\ of the numerically obtained 
transmission coefficient T num as a function of Kd. Here, TV = 1000 
segments are used. The corresponding barrier is shown in Fig. |2la), the 
effective mass is assumed to be constant. 



coefficient T, giving the tunneling probability of an electron 
[4 |. The unsymmetrized, symmetrized and Airy function trans- 
fer matrix approaches are evaluated, based on the expressions 
dD, (fl2] | and dlQ) . respectively; for comparison, also the 
QTBM result is computed. Assuming an electron energy of 
E = and a constant effective mass of m* = 0.067 m e 
corresponding to GaAs, where m e is the electron mass, the 
exact value obtained by evaluating (fl4l > is T = 1.749 x 10~ 4 . 
Fig. |3] shows the relative error e T (N) = |1 - T num (A) /T\ 
as a function of A oc A" 1 . Here, T num (iV) is the numerical 
result for the transmission coefficient, as obtained by the 
different methods for a subdivision of the structure into N 
segments of equal length A = d/N oc A^ 1 . As can be 
seen from Fig. [3] the error scales with A^ 1 cx A for the 
unsymmetrized transfer matrix approach and with A~ 2 oc A 2 
for the other methods. This can easily be understood by means 
of the local discretization error, which is O (A 3 ) for the Airy 
function approach and the symmetrized transfer matrix, and 
O (A 2 ) for the unsymmetrized matrix, as discussed above 
and in the appendix. When the overall transfer matrix of the 
structure is computed from (J8]l, the individual LDEs arising for 
each of the N segments accumulate, thus resulting in a total 
error NO (A 2 ) = O (A) for the unsymmetrized approach and 
NO (A 3 ) = O (A 2 ) for the other schemes. 

The symmetrized transfer matrix approach and the Airy 
function method are the most accurate, both exhibiting a 
comparable error et (A). However, the symmetrized matrix 
approach is much more computationally efficient, being over 
20 times faster than the Airy function method in our MATLAB 
implementation. For a given A, the QTBM is even three times 
faster than the symmetrized matrix approach, but also 40 times 
less accurate, meaning that it requires \/40 s» 6 times as many 
grid points as the symmetrized matrix approach to achieve the 
same accuracy. 

Fig. shows again the relative error et, but now for 
a fixed number of segments A = 1000. Instead, the shape 
of the potential is modified by varying K in (fl3l >, and also 
adapting Vb and V\ so that V (z) remains constant at z = 
and z — d and only the curvature of the potential changes. 




12 3 4 

10 10 10 10 

Number of segments N 

Fig. 5. Relative error ee = 1 — Enum/E\ of the numerically obtained 
eigenenergy E nnia for the (a) first and (b) second bound state as a function 
of the number of segments TV. The corresponding well is shown in Fig.ffjb), 
the effective mass is assumed to be constant. 



The symmetrized matrix approach and the Airy function 
method exhibit a superior accuracy especially for small K, 
corresponding to a weak curvature of the potential. While the 
error of the unsymmetrized matrix approach and the QTBM 
show only a weak dependence on K, the Airy function method 
becomes exact for K — > 0, where the potential becomes 
piecewise linear. Interestingly, also the symmetrized transfer 
matrix approach has a vanishing error et for a specific value 
of K, at Kd& 0.167. 

Now we apply the different numerical methods to the bound 
states of the potential well shown in Fig.[2fb). Again assuming 
a constant effective mass of m* = 0.067 m e , evaluation of (fl4l 
yields two bound states with eigenvalues E\ = — 0.1343eV 



ACCURACY OF TRANSFER MATRIX APPROACHES FOR SOLVING THE EFFECTIVE MASS SCHRODINGER EQUATION 



5 



and Ei — — 0.0129 eV, respectively. In the following, we 
compare the accuracy of the numerically found eigenenergies 
-E-num, as obtained by the unsymmetrized and the symmetrized 
transfer matrix approach and the Airy function method, cor- 
responding to the expressions ©, ( fT2l and JTOt . respectively. 
Here, we again divide the structure into N segments of equal 
length A = d/N oc TV -1 . Fig. [5] shows the relative error 
£e{N) = |1 - E num (N) /E\ for the first and the second 
bound state as a function of AT. As for the transmission 
coefficient in Fig. [3] the error scales with iV -1 oc A for 



the unsymmetrized matrix approach and with N 



the other methods. Again, the symmetrized matrix approach 
and the Airy function method exhibit a comparable value 
of £e (-/V), being far superior to the unsymmetrized matrix 
approach. 

B. Position dependent effective mass 

Now we compare the accuracy of the different methods for 
a position dependent effective mass m*(z). Here, we choose 
the same exponential ansatz for the potential as above, see 
( fT3] l and Fig. [2] For an effective mass of the form m* = 
uiq exp (— Kz), again an analytical solution exists: 

ip = ci J M (a) exp (-Kz/2) + c 2 Y„ (a) exp {-Kz/2) , (16) 

with 



a(z) = 2 



1 



, TO Vi 



^-2m (V - E) 
HK 



exp 



— Kz 
2 



(17) 



The transfer matrix definitions (O, (TTZt are also valid for 
position dependent effective masses. In the Airy function 
approach ( [Tol l, a position dependent effective mass can be 
accounted for by assuming a constant value within each 
segment, as discussed at the end of Section IH-AI Here, we 
assign the averaged mass (m* +m* +1 ) /2 rather than m* 
to each segment, since then the third order LDE, found for 
the amplitudes A and B in the case of position independent 
masses, is also preserved for the position dependent case, see 
the appendix. 

Fig. [6] corresponds to Fig. [3] but now for a position 
dependent effective potential with m* — 0.2 m e exp (— Kz) 
for < z < d and m* = 0.067 m e otherwise. The exact 
transmission coefficient for an electron with energy E = 0, 
as obtained by evaluating ( fl6b . is now T = 5.376 x 10~ 10 . 
From Fig. [6] we can see that also here the error scales with 
AT -1 ex A for the unsymmetrized matrix approach and with 
A^~ 2 ex A 2 for the other methods, compare Fig. [3] Again, the 
symmetrized matrix approach and the Airy function method 
are the most accurate, with the symmetrized matrix approach 
being numerically much more efficient. 

For the sake of completeness, Fig. (Q is shown as the 
counterpart of Fig. 2J but now taking into account a position 
dependent effective mass as above. Again, the symmetrized 
matrix approach and the Airy function method have a superior 
accuracy especially for small values of K, corresponding to a 
weak curvature of the potential. 
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The corresponding barrier is shown in Fig. 12 a), the effective mass is assumed 
to be position dependent. 
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Fig. 7. Relative error et = |1 — T nuul /T\ of the numerically obtained 
transmission coefficient T nU m as a function of Kd. Here, TV = 1000 
segments are used. The corresponding barrier is shown in Fig. [2ja), the 
effective mass is assumed to be position dependent. 



IV. Example: Schrodinger-Poisson solver 

In the following, we apply the transfer matrices discussed 
above to a real-world example, namely finding the wavefunc- 
tions and eigenenergies of the quantum cascade laser (QCL) 
structure described in |fT31 . The goal is to evaluate and com- 
pare the performance of the different approaches for a practical 
problem, and to discuss the inclusion of additional important 
effects. Specifically, we here also account for energy-band 
nonparabolicity, and complement the Schrodinger equation by 
the Poisson equation to take into account space charge effects. 
In practice, extensive parameter scans have to be performed 
for QCL design optimization. Thus, the simulation of QCLs 
calls for especially efficient methods, the more so as the self- 
consistent solution of the Schrodinger-Poisson system results 
in a further increase of the numerical effort. 

In simulations, the QCL structure is defined by an infinitely 
repeated elementary sequence of multiple wells and barriers 
(called a period). For such a structure under bias, it is 
sufficient to compute the eigenenergies and corresponding 
wave functions for a single energy interval given by the bias 
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across one period; the solutions of the other periods are then 
obtained by appropriate shifts in position and energy. We 
solve the Schrodinger equation using the approaches defined 
by ©, ( TTZi i and ( [Tol l, respectively. For all three methods, we 
treat band edge discontinuities at the barrier-well interfaces 
explicitly using the matching conditions (0), to obtain an 
optimum accuracy. We use a simulation window of four 
periods to keep the influence of the boundaries negligible, 
and determine the bound states similarly as in Section [HI] To 
combine reasonable numerical efficiency with a good accuracy, 
we choose a segment length of A = 2nm (the last segment 
of each well or barrier is A < 2 nm). 

Various models are available for including nonparabolic- 
ity lfl6l : here, we use an energy dependent effective mass 
to* (E) = m* [1 + (E- Vj) /E g>j ] (with band gap energy 
E g j at position Zj), which can straightforwardly be imple- 
mented into the transfer matrices. The Poisson equation is 
given by 0, El 



Unsymmetrized matrix (2 nm) 
Symmetrized matrix (2 nm) 
Airy functions (2 nm) 



d z [e (z)d z <p(z)] = e 



(18) 

leading to an additional potential -e<p in ([T). Here, e (z) is the 
permittivity, e is the elementary charge, N (z) is the doping 
concentration, and ri2D,n is the electron sheet density of level 
n with wave function ip n (z). While for an operating QCL, 
«2D,n can only be exactly determined by detailed carrier 
transport simulations Q, this is prohibitive for design opti- 
mizations of experimental QCL structures over an extended 
parameter range. Thus, for solving the Schrodinger-Poisson 
system, simpler and much faster approaches are commonly 
adopted, such as applying Fermi-Dirac statistics 0, ifTTIl 



«2D,n 



fceTln ^1 + exp 



H - E n / (k B T) 



(19) 



where /i is the chemical potential, fee is the Boltzmann 
constant, T is the lattice temperature, and to* is the effective 
mass, here taken to be the value of the well material. In 
( fT9l . we use the energy of a state relative to the conduction 
band edge E n — E n — J V |"0 n | 2 dz rather than E n itself 
to correctly reflect the invariance properties of the biased 
structure. Especially, this ensures that the simulation results 
do not depend on the choice of the elementary period in the 
structure. The chemical potential /j, is found from the charge 
neutrality condition within one period. The Schrodinger and 
Poisson equations are iteratively solved until self-consistency 
is achieved. For the Poisson equation ( fT8l , we employ a finite 
difference scheme on a 1 A-grid, where we use (fZ|i and (O to 
appropriately interpolate the eigenfunctions obtained from the 
Schrodinger solver. 

Simulations of the QCL in [15] have been performed at 
various temperatures, considering the seven lowest levels (i.e, 
with lowest energies E n ) within each period. In Fig. ||8), 
the obtained energy levels and wave functions squared of a 
single period are shown for the unsymmetrized, symmetrized 
and Airy function matrix approach at T = 300 K, using a 
segment length of A = 2nm. Also the symmetrized transfer 
matrix result for A = 0.1 nm is plotted for reference. The 




Fig. 8. Self-consistent band profile (grey line), energy levels and wave 
functions squared for the QCL in 1151 at a bias of 48 kV /cm and a 
temperature of 300 K. Shown are the results as obtained with the three transfer 
matrix methods for A = 2 nm, and also the symmetrized matrix result for 
A = 0.1 nm, which practically coincides with the symmetrized matrix and 
Airy function results for A = 2 nm. 



symmetrized matrix and Airy function results exhibit a similar 
accuracy, with deviations in eigenenergies of around O.lmeV 
from the high-accuracy result obtained with A = 0.1 nm. In 
Fig. ©, those three curves are practically indistinguishable. 
On the other hand, the unsymmetrized matrix method pro- 
duces deviations of around 5meV. The unsymmetrized and 
symmetrized matrix approach require approximately the same 
computation time for obtaining the self-consistent result in 
Fig. (JHJl, while the Schrodinger-Poisson solver based on the 
Airy functions is about 10 times slower. This confirms that the 
symmetrized transfer matrix method combines high numerical 
efficiency with excellent accuracy for practical applications. 

V. Conclusion 

In conclusion, we have compared the accuracy of different 
transfer matrix approaches, as used for solving the effective 
mass Schrodinger equation with an arbitrary one-dimensional 
potential and a constant or position dependent effective mass. 
In particular, the local discretization error has been derived for 
the Airy function approach resulting from a piecewise linear 
approximation of the potential, and for unsymmetrized and 
symmetrized transfer matrices based on a piecewise constant 
potential approximation. Furthermore, numerical simulations 
have been performed to evaluate the numerical accuracy of 
the different approaches for scattering and bound states, em- 
ploying exponential test potentials. Comparisons to the finite 
difference method, specifically the QTBM, have also been 
carried out. Additionally, self-consistent Schrodinger-Poisson 
device simulations are presented. 

The symmetrized transfer matrix approach and the Airy 
function method exhibit a comparable accuracy, being superior 
to the other methods investigated. However, the symmetrized 
matrix approach achieves this at a significantly reduced nu- 
merical cost, moreover avoiding the numerical problems asso- 
ciated with Airy functions. All in all, the symmetrized transfer 
matrix approach is shown to combine the numerical efficiency 
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and straightforwardness of its unsymmetrized counterpart with 
the superior accuracy of the Airy function method. 

Appendix A 
Local discretization error 

In the following, we derive the local discretization error 
(LDE) for the different types of transfer matrices. In this 
context, we investigate the piecewise constant potential ap- 
proximation based on matrix (f5]l and its symmetrized version 
( fl2] >. as well as the piecewise linear potential scheme dTOb . 
As mentioned in Section |IlJ the segments are chosen so that 
band edge discontinuities in the structure coincide with the 
borders between two segments, enabling an exact treatment 
in terms of the matching conditions ([3]). Thus, for our error 
analysis we imply that the potential and effective mass vary 
smoothly within each segment, i.e., have a sufficient degree 
of differentiability. Otherwise, no further assumptions about 
the potential shape and effective mass are made. The local 
discretization error for tp at Z — Zj^\ IS 

rf +1 = tpj+x - ip Oj+i) , (20) 

where ipj+i is the approximate wavefunction value at Zj+i 
obtained by the transfer matrix approach from a given value 
tp (zj) — ipj at Zj, while tp (zj+i) is the exact solution. For 
evaluating the LDE, it is helpful to express ip (zj+i) in terms 
of a Taylor series, 

tP (z j+1 ) = ^ + tp'jAj + i$'A» + gVf A? 

+ ^f A i+<?( A i)- (21) 
Analogously, we can define an LDE for the derivative ip', 

rf +1 = $+i - (z j+ i) , (22) 
and express ip' (zj+i) as 

= iP' j + iP'jAj + \tpf ] A] + AJ + O (A|) . 

(23) 

A. Piecewise constant potential approximation 

Using (O, we can relate Aj and Bj to the wavefunction at 
position zj, 

^=K^ +i ^0' (24) 

and express the LDEs for the amplitudes A and B in terms 

of t/ +1 and rf +1 , 

T f+i = A 3+l - A ( z j+l) = \ (rf+i - i fc^ T J+ 1 ) ' 
tf+i = B J+ i B (z J+1 ) = \ [rf +l + i^^i) • (25) 



For the unsymmetrized transfer matrix, we obtain from (f2| 
with the expressions and Q 

tpj+i = Aj + i + Bj +1 = exp (ikjAj) Aj + exp (— ikjAj) Bj, 
tp' j+1 = ifc i+ i [A j+1 - B j+1 ) 

= ifcj+i ^ J [Aj exp (ikjAj) — Bj exp (— ikjAj)] . 
Pj+i 

(26) 

For calculating the LDE, we insert the expressions d26| i and 
( ETT i into (|20l l, where we express Aj and Bj by (124-b and use 
(Q~|l to rewrite the derivatives ipj in ( ETT i as 

-Ml 

t ' =- k ^j + ^'j, 
J J J m * J 

4 ) =-^^-^' j + ^' j , (27) 
with kj = ^J2nij (E — Vj ) jh. A Taylor expansion then yields 

T f+i = cos (kjAj) ipj + sin (kjAj) kj l tp'j 

-^-^A'-^Al-^fAf + OC^) 

lmf . - if (m*B)' to*" \ 
2 TTij 6 \ to* m* J I J 

+ O (Af) . (28) 

Analogously, by inserting the expressions d26l i and d23l into 
d22T i we obtain 

3 

-rP'j-^Aj-^fA" j+ 0(A^ 

( 771* j_i 777*' \ , 777* — 777* , , 

\ TO* TO* / J 777* 

1 / (m*k^)' TO*" , 777* — TO* , -i \ 

2 \ 777* 777* 7 TO* y J 

+ 0(AJ) 

-(^-^)* A 5 + °( A 5)- < 29 » 

where we use to* +1 = to* + mfAj + \m*" A 2 j + O (A|) 
to obtain the last line of ©. Thus, r/ +1 is O (Af) (O (A|) 
for a constant effective mass), and Tj +1 = O (Aj). With ( f25l l, 
we see that both 7^j. x and r£. x are O (A|). 

In a similar manner, we obtain for the symmetrized matrix 
^ T f +l = O (Af) and rf^ = O (A?). More precisely, the 
calculation yields for a constant to* 

<i - ( fc f)"^ A ? " ^ (^)'' '^ + ° ( A D 00) 

(and a somewhat more complicated expression for a position 
dependent effective mass). This means that tA_ 1 and t? +1 are 
now O (Af). 
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B. Piecewise linear potential approximation 

For computing the LDEs ( f20b and (1221 of the Airy function 
approach, we proceed in a manner similar as above. Equations 
( [Tol l and dTTb give the relation between values ipj, i// at Zj 
and the numerical result ipj+i, V'j'+i obtained at Zj+i from 
the Airy function approach: 



1 



D , 



Ai( Sj + ^)Bi'( Sj ) - Ai'( Sj )Bi( Sj + ^) 
Ai( flJ )Bifo + -/) - Ai( Sj + ^)Bi( Si ) 



Is, 



1A? 



1 s 



^ + i|A3 + l|)^ + (A^) 



A, 



(31) 



Ai'( Sj + -/)Bi'( Si ) 



-Bi'fo + -/)Ai'( Sj -) 



AifoJBi'fo + ^) - Ai'fo + ^)Bi( Si ) 



Ai 1 1 



is 



'i + i^ + 5f)« + o(Aj), 



(32) 



with = Ai(sj)Bi'(sj) — Ai'(sj)Bi(sj). The exact results 
(Zj+i) and if)' (zj+%) are again expressed by the Taylor 
series expansions (l2TT i and d23l . respectively, where we rewrite 
the derivatives tpj in terms of ipj and ip'-. For a constant 
effective mass, we have 



ib (3 > = £T 2 Si ip' 4 + £- 3 ^-ib 



(4) _ ,-4 2, 



VI' 



3 r J 



v z . 



(33) 



with V^.j = {Vj+i — Vj) I Aj, and obtain with the expressions 
©, (ED and (EB 

1 



i = ^+1-^(^+1) 



j + L Tj-r± r v-jTi, 24 

and with <|22), (|23) and d32j 



(fcf ) ^ A|+0 (A|) , (34) 



rf +l = - ^ fe+i) = (fc, 2 )" ^A" + ° (A) . 

(35) 



where kj = J 2m* (E — Vj) / K. Using ( fTTT i. we can express 
the LDEs for the amplitudes A and B in terms of Tj' +1 and 

rf +1 , and obtain rf +1 = O ( A?) rf +1 = O (Af) . 

As described in Section III-A1 a position dependent ef- 
fective mass can in the Airy function approach be treated 
by assuming a constant value within each segment j, e.g., 



m* = to* (Zj), and applying the matching conditions <[3j at 
the section boundaries [2|. The result for t/'j+i m d32l has 
thus to be multiplied by to* +1 /to* before inserting it into 
03. While r/^ is still (A|), t/ +1 drops to O (Af), now 
yielding rf +1 = O (A|), rf +1 = O (A|). The error analysis 
also shows that and thus T j+i, T ^+i can be improved 

to O (A|) by assigning an averaged mass (to* + to* +1 ) /2 
rather than to* to each segment, and applying the matching 
conditions correspondingly. 
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